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ABSTRACT 

Density waves excited by planets embedded in protoplanetary disks play a 
central role in planetary migration and gap opening processes. We carry out 2D 
shearing sheet simulations to study the linear regime of wave evolution with the 
grid-based code Athena, and provide detailed comparisons with the theoretical 
predictions. Low mass planets (down to ~ 0.03M e at 1 AU) and high spatial 
resolution (256 grid points per scale height) are chosen to mitigate the effects 
of wave nonlinearity. To complement the existing numerical studies, we focus 
on the primary physical variables such as the spatial profile of the wave, torque 
density, and the angular momentum flux carried by the wave, instead of secondary 
quantities such as the planetary migration rate. Our results show percent level 
agreement with theory in both physical and Fourier space. New phenomena such 
as the change of the toque density sign far from the planet are discovered and 
discussed. Also, we explore the effect of the numerical algorithms, and find that 
a high order of accuracy, high resolution, and an accurate planetary potential are 
crucial to achieve good agreement with the theory. We find that the use of a too 
large time-step without properly resolving the dynamical time scale around the 
planet produces incorrect results, and may lead to spurious gap opening. Global 
simulations of planet migration and gap opening violating this requirement may 
be affected by spurious effects resulting in e.g. the incorrect planetary migration 
rate and gap opening mass. 

Subject headings: Planet-disk interactions, Protoplanetary disks, Hydrodynam- 
ics, Methods: numerical, Planets and satellites: formation 
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INTRODUCTION 



The discovery of a number of "hot Jupiters" - - extrasolar giant planets residing on 
very tight orbits — calls for understanding the dynamical pathways that brought these ob- 
jects to their current orbits. Due to the difficulty of forming these planets in situ, it was 
suggested that potential hot Jupiter's form at larger separation s from their stars and are sub- 



sequently transported inwards by some migration mechanism (jPapaloizou fc Terqueml 12006 



Papaloizou et al.l 120071). In protoplanetary disks, p lanets act as sources of nonaxisymmetric 
density waves ( IGoldreich &: Tremaindll979l Il980l hereafter GT80), and coupling of plane- 
tary gravity to these waves acts to transfer angular momentum between the disk and the 
planet. This ultimately results in the orbital migration of embedded planets. At the same 
time, density waves generated by planets can affect the structure of the disk (potentially 
leading to gap opening) once the angular momentum that they carry is transferred to the 
disk material. Thus, better understanding of t he density wave properties and processes of 
their excitation, propagation, and damping plays a crucial role in the global picture of planet 
formation. 

Evolution of the density wave can be separated into the linear and nonlinear stages. The 
former applies to density waves excited by low-mass planets, which have not had enough time 
to propagate very far from their launching site (a more precise condition is formulated in jQ. 



The linear theory of disk-planet interaction was pioneered by IGoldreich fc Tremaind ( 11979 



GT80) who provided a description of the density wave excitation in two-dimensional (2D) 
gaseous disks. This theory has been subsequently refined in various ways, e.g. by calculating 
the asymmetry of planetary interaction with inner and outer portions of the disk (Ward 
1986), and extending it to fully three-dimensional disks (ITanaka et al.ll2002l ). 



These studies concentrated on understanding th e behavior of perturbed quanti ties Fourier 
decomposed in azimuthal angle ( )Artymowiczlll993al Jbl; lKorycansky fc Pollacklll993l ). often ad- 
dressing only the integral wave properties such as the total angular momentum carried by 
the waves, or the differential Lindblad torque (jWardl Il986l ). Only recently has the linear 
evolution of the wave properties in physical space, namely th e spatial dependence of the 



perturbed fluid velocity and surface density, been addressed bv lGoodman &: Rafikovl (12 001 



hereaf ter GR01) in the loc al shearing sheet approximation. Subsequently lOgilvie fc Lubow 
( 2002 ) and Rafikov (j2002a ) studied the shape of the density wake in the global s etting. This 



paved way for subsequent exploration of the nonlinear wave evolution (GR01, iBate et al. 
20021 ) . understanding of which requires knowledge of the wave characteristics in physical 
space. 



In parallel with these analytical and semi-analytical studies, numerous attempts to ver- 
ify linear theory by direct hydrodynamical simulations were made. However, the majority 
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of these studies compared with theoretical predictions only the secondary or derived prop- 
erties of the planet-disk syst em related to the densit y wave excita t ion and evolution. Fo r 
exa mple, in their simulations ILin fc Papaloizoul (Il993h , IWardl ( 119971 ) , ID'Angelo et al.l ( 120031 ) , 
and ID'Angelo fc Lubowl (120081 ) measure the rate of planet migration due to its gravitational 
interaction with the disk — a quantity that is determined by the small difference between 
the summed (over all azimuthal harmonics) one-sided t orques that the planet exerts on the 



inner and outer portions of the disk. Other studies ( ILin fc Papaloizoul Il993t iMuto et al. 
20101 ) tried to measure the planetary mass at which a gap opens in the disk — a process that 



depends not only on the magnitude of the angular momentum carried by the waves but also 
on the intricate details of their dissipation. In both situations the verification of the linear 
theory ends up being rather indirect. 

It is certainly much more convincing and more robust to verify linear theory by compar- 
ing with its predictions the primary wave properties, such as the spatial distribution of the 
perturbed fluid properties deri ved from simulations . Rece ntly some effort in this direction 
has been made. In particular, ID'Angelo fc Lubowl ( 120081 ) investigated the dis tribution of 



the to rque imparted by the planet on the disk as a function of radial distance. IMuto et al. 



( 120101 ) showed the spatial structure of the density perturbation associated with the wave 
obtained in their 2D simulations. Unfortunately, in both cases no quantitative comparison 
with linear theory predictions was provided. 

In this work, we use 2D local inviscid hydrodynamical simulations in the shearing sheet 
geometry to verify the linear theory. For the first time we provide direct comparison of 
the primary physical variables — perturbed surface density profile, torque density in both 
physical and Fourier space, and angular momentum flux — derived from simulations with 
the theoretical predictions. Our verification effort is the most thorough, sophisticated and 
quantitative to date, and we achieve agreement with theory at the level of several per cent 
for all physical variables explored. In addition, using our numerical results we address some 
remaining issues in the linear theory , such a s the different expressions for the torque density 
in Fourier space (GR80 jArtymowiczlll993al ). We also formulate a set of requirements which 
should be satisfied by simulations to reach convergence and to reliably reproduce the linear 
theory results. The way in which we make the comparison an d the agreement we achieve 
can serve as a standard for checking the performance of different codes. 

The paper is structured as follows. We briefly summarize the main results of the linear 
density wave theory in £j2j A description of the code and our numerical setup are given 
in §21 including a discussion of potential problems associated with implementing orbital 
advection algorithm in studies of the disk-satellite interaction. We present our numerical 
results and compare them with theoretical predictions in §U followed by the investigation of 
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the sensitivity of simulation outcomes to the different numerical algorithms in £j5j In §|6] we 
describe and discuss a commonly ignored numerical problem in disk simulations which if not 
handled appropriately can lead to incorrect results. Finally, a summary is provided in §7J 



2. Summary of linear theory results. 



In linear theory of disk-planet interaction laid out in cylindrical geometry (IGoldreich fc Tremaine 



19791 . GT80), the planetary potential is normally decomposed into a series of Fourier har- 
monics. Assuming all perturbed fluid variables to be proportional to exp[i(m(ip — Q p t) + 
J r k r (r')dr')] (m is a positive integer) and neglecting the self-gravity of the disk, one can 
obtain the dispersion relation for the m-th wave harmonic in the WKB limit in the form 



m 



Mr) - n p ] 



2^2 



+ k;c 



(1) 



where k is the epicyclic frequency (k, = Q for a Newtonian potential), k r is the wavenumber, 
Q is the fluid angular frequency, and Q p is the pattern speed of the perturbation equal to the 
angular frequency of the point-mass perturber moving on a circular orbit with semi- major 
axis r p , that is Q p = tt(r p ). As Eq. flTJ demonstrates, there are two (inner and outer) 
Lindblad resonances (LR) where k r — » and m-th harmonic of planetary potential couples 
best to the wave-like perturbation. These resonances are located at radii ri,rn given by 



Q p m/(m ± 1), which in the limit of m ^> 1 results in 



XL,: 



±- 



2h 
3/i' 



(2) 



where x = r — r p , h 



c s /Q p is the vertical disk scale length, and fi = mh/r p . 



The superposition of all the harmonics with different azimuthal wavenumbers m forms 
inner and outer spiral density waves. In physical space fluid perturbation associated with 
the wave is concentrated in narrow wake (with azimuthal extent of order h), the location of 
which is described in the local limit by (GR01) 



Vv 



— sgn(x^ 



3^ 
4 h 



(3) 



where y = r p {ip — <p p ) and azimuthal angle (p is counted in the prograde sense (with respect 
to planetary orbital motion). This location corresponds to the stationary phase of different 
Fourier harmonics composing the wave, as can be easily seen from Eq. ((T|). The wake shape 
in the global setting accounting for cylindrical geo metry of t he disk and Keplerian profile of 
n(r) was explored bv lOgilvie fc Lubowl ( 120021 ) and iRafikovl (j2002al ). 
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Each of the perturbation harmonics carries angular momentum flux, which for the m-th 
harmonic is given by (GT80). 



F H (m) = F^(m)^^- y (4) 
F WKB (m) = 4m^o / GM,V m + )]3 _ (5) 



3 ^ V r p 7 

Here S is the unperturbed surface density of the disk, F^ KB {m) is the asymptotic Fourier 
contribution of the flux reached beyond the launching region on one side of the disk calculated 
in the WKB approximation (K and K\ are the modified Bessel functions of zeroth and first 
order), and the expression is only accurate for m 3> 1. The ratio FH{m)/F^ KB (m) accounts 
for the difference between the precise value of Fn(rn) and F^ KB (m). This difference becomes 
important for m > r p /h 3> 1 since, as demonstrated by GT80, Fh(iti) decays exponentially 
with m at high azimuthal wavenumbers — the so-called "torque cutoff" phenomenon. As a 
result, the total angular momentum flux carried by the wave 

Fh = Y, Mm) « 0.93 (GM P ) 2 ^f^, (6) 

m=l s 

is dominated by harmonics with m ~ r p /h (or /i ~ 1), which are excited within a radial 
separation of order h from the planet, see Eq. $2§. This local nature of the wave excitation 
has important implications for the possibility of analytical treatment of the wave evolution. 

In real protoplanetary disks, the torques acting on the inner and outer disks are slightly 
different (fractional difference ~ h/r p <gj 1). givin g rise to a net torque on the planet and 



leading to Type I planetary migration (jWardll 19861 ). However, in shearing sheet configuration 



the inner and outer parts of the disk are identical so there is no net torque acting on the 
planet and migration does not arise. 

There is an important planetary mass scale characterizing excitation and propagation 
of the density waves, the so called thermal mass (GR01), which we defkjj] here as 

There are several reasons for the significance of this mass scale. First, a Hill radius of the 
planet with M p = M t j, is of order the scale height of the gaseous disk h. Second, the Bondi 
radius Rb = GM p /c 2 s , which is the size of a region in which planetary gravity strongly 



lr This definition differs from that of GR01 who used Mi = 2cl/{Z%G) = (2/3)M th instead of M, 



th- 
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perturbs local pressure distribution, also becomes comparable to h for M p = M t h (Rafikov 
2006). Third, and most importantly for our present study, planets with M p > M th generate 
density waves which are nonlinear from the very start (<5£/£ ^ !)• 

In the small planet mass limit of 



perturbations induced by the planetary gravity are weak (5S/E < 1) and excitation of the 
density wave by the planet can be studied purely in the linear approximation. It is in this 
limit that the linear theory of wave driving developed by GT80 and all of its results are valid. 
In the opposite limit of M p > M t h the nonlinear effects must be accounted for even at the 
stage of wave excitation, which is not feasible analytically. 

As mentioned above, the generation of the wave is largely completed at a distance of 
order h away from the planet, because this is where the most important Lindblad resonances 
lie. After that, in the absence of viscosity, the wave travels freely, no longer affected by 
the planetary gravity. At this stage of propagation nonlinear effects start to accumulate and 
finally result in wave breaking, formation of a shock, and transfer of wave energy and angular 
momentum to the disk fluid. Prior to the steepening into a shock the angular momentum 
flux Fh carried by the wave far from the planet (outside of the excitation region \x\ ~ h) is 
exactly conserved in the absence of linear dissipation. This nonlinear stage has been explored 
analytically in GR01 in the limit OH]). The numerical results on the nonlinear wave evolution 
will be presented in ?, here after Paper II. 

We do not cover this stage in this work, but we mention the result of GR01 for the 
radial separation away from the planetary orbit at which the shock induced by the nonlinear 
effects first occurs: 



where 7 is the adiabatic index of gas. For M p <C M th one finds the shock to be well separated 
from the planet, Z s h > h, so that the wave excitation occurs in the linear regime, as asserted 
before. However, even for M p -C M t h the linear theory of wave propagation fails for x > l s ^. 
In this work we are always exploring the low planetary mass limit (jSJ) and consider only 
I a; I < l s h so that linear theory applies. 



M p <M t 



th 



8 




(9) 



3. Numerical setup 



Most analytical work on the linear evolution of density waves was carried out in ap- 
plication to two-dimensional (2D) disks. To provide a meaningful comparison with these 
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studies, in our current work we employ 2D hydrodynamical simulations in local shearing box 
configuration. The computational tool for this investigation is Athena, a grid-based code 
for astrophysical gas dynamics using higher-order Godunov methods. Athena is written in 
conservative form, so it conserves mass, momentum, and energy d own to machine precision 



The mathematical foundations of the algorithms are described in iGardiner fc Stone! ( 12005 



20081 ). a nd a comprehensiy e description of the implementation and tests of the algorithms is 



given in IStone et al.l (120081 ) . 



The implementation of the shearing box approximation in Athena is described in lStone fc Gardiner 



( 120 lOh . This approximation adopts a frame of reference located at radius r p corotating with 
the disk at orbital frequency Q p = Q(r p ). In this frame, the 2D hydrodynamics equations 
are written in a Cartesi an coordinate system (x, y, z) that has unit vectors i, j, and k as 
rtstone fc Gardiner! boioh : 



f + V.(Ev) 
<9£v 



0. 



+ V-(Svv +p) = m 2 p (2qxi) - 2Q p k xEv- SV$ P , 



dt 
OE 

— + V ■ (Ev + p ■ v) = fiJSv ■ (2qxi) - Sv ■ V$ p , 



(10) 

(11) 
(12) 



where $ p is the gravitational potential of the planet (whose form we will discuss in the next 
section), p is the gas pressure, E is the total energy density (sum of the internal and kinetic 
energy, but excluding the gravitational energy). The shear parameter q is defined as: 



dlnr 



(13) 



so that for Keplerian flow q 
is the ideal gas law with 



3/2. The equation of state (EOS) we use in the simulations 

£in = ~~~T ( 14 ) 

7-1 

where E- m is the internal energy and 7 = 5/3. We also use an isothermal EOS, in which 
case p = SCg, where c s is the isothermal sound speed. We do not include explicit viscosity 
(i.e. there is no explicit linear dissipation in the system) and we do not account for the 
self-gravity of the disk (thus we concentrate on studying the low mass disks). 

To accurately resolve the flow in the vicinity of a gravitating mass without making the 
time step too short one usually resorts to softening the potential of the perturber. There 
are many ways of doing this but all of them come up with the form of the potential that 
converges to the Newtonian potential = —GM p /p at large separations p = \J x 2 + y 2 , 
beyond the softening length r s . At small separations, p r s smoothed potentials become 
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finite allowing the time step of simulations to stay finite. To test the sensitivity of our results 
to the specific method of potential softening we tried three different forms of the softened 
planetary potential. The second order potential 

^ = -^P (p2 + ,2)1/2 ( 15 ) 

converges to at p ^> r s as (r s / p) 2 (which means the fractional error is 0((r s /p) 2 ) as 
r s /p — > 0). This is a standard form of the potential used in a majority of numerical hydro- 
dynamic studies. We also studied the fourth order potential 

converging to the point mass potential as (r s /p) for p 3> r B , and the sixth order potential 

$( a) = _ GM P 4 + 5 As72 + 15r a V8 

(p2 +r 2)5/2 ^'J 

converging to $k at p ^> r s as (r s / p) 6 . The difference in accuracy with which these potentials 
represent &k at large p can be quantified by the distance from the planet at which a given 
potential deviates from $^ by 1%. This distance is 7r s for $ p 2 \ 2.3r s for $p 4 \ and 1.5r s for 
$ p 6 ' ) . Note that the first-order potential = — GM p /\p + r s \, which we do not use here, 
gives only 10% accuracy relative to $x even at p = 10r s , and we strongly discourage its use. 
The effects of different potential prescriptions will be discussed in detail in §5.31 



To guarantee accuracy and convergence of our results we carry out an exhaustive ex- 
ploration of various numerical parameters characterizing our simulations, described in §5j 
Unless noted otherwise, our results shown in §H are obtained with the following numerical 
parameters: an isothermal equation of state, the Roe solver with third order reconstruc- 
tion in characteris tic variables, and the corner transport upwind (CTU) unsplit integrator 



(IStone et al.l 120081 ) , a resolution of 256 grid points per scale length h (subsequently denoted 
256//z for brevity), and fourth order accurate potential $ p 4 ^ given by Eq. ( flBT) with softening 
length r s = h/32. 

High order of accuracy, high resolution and an accurate form of planetary potential used 
in our runs allow us to properly capture the details of wave evolution. High resolution ensures 
low levels of numerical viscosity and prevents the angular momentum accumulated by the 
wave from spurious dissipation (see Section [5] for more discussion). We run a series of test 
simulations with otherwise identical conditions but different explicit Navier-Stokes viscosity, 
and measure the wave properties. As the explicit viscosity decreases the simulation results 
gradually converge to the one with zero explicit viscosity, which indicates the numerical 
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viscosity dominates the explicit viscosity. For a typical simulation with low M p = 2.09 x 
10~ 2 M t h and an isothermal equation of state, the effective Shakura-Sunyaev a-parameter 
(a = u/hcg) characterizing our numerical viscosity is found to be b elow 10~ 5 - Suc h small 
levels of viscosity are expected in dead zones of protoplanetary dis ks (Gammiel 199~ti), wher e 



magnetorotational instability (MRI) may not operate effectively ([Fleming fc Stone! 120031 ). 
Also note even when the MRI operates, MHD turbulence may not act like a Navier-Stokes 
viscosity. 

For the standard simulations in this work to fit the parabolic wake (see Eq. [3]) in the box 
we use box size 12h x 6 Ah, thus the overall grid resolution in our runs is about 3072 x 16384 
(for the standard resolution of 256/ h). In a few cases with very small M p , we extend the 
simulation box size to 20h x 156/z to trace linear wave evolution further out in x, since smaller 
M p delays shock formation (see eq. OH])). Our simulations are run for at least 10 and in some 
cases up to 50 orbital periods. Figured] presents a typical snapshot of one of our simulations 
showing the density structure and the spiral waves. 

We use the following boundary conditions (BCs). On x (radial) boundaries, we keep 
values of all physical variables in ghost zones fixed at their respective unperturbed Keplerian 
values {i.e. keep the ghost zones as their initial states), and the waves leave through the x 
boundaries when they reach the edge. In our shearing sheet simulations, experiments show 
that this x BC has less wave reflection and outflowing of the fluid than the conventional 
outflow x BC, in which case the ghost zones are copied fr om the last activ ely-updated 



column of cells of the simulations at every time-step (see also lMuto et al.ll2010l ). We do not 
expect our adopted radial BC to affect wave evolution since significant radial fluid motions 
are not expected to arise in our simulations anyway. 

On the y (azimuthal) boundary, we experimented with two BCs: the conventional out- 
flow BC, as described above, and an inflow/outflow BC. In the latter case, the variables in 
the ghost zones are fixed at their initial values if they are the physical "inflow" boundaries 
(the regimes y < 0,x < and y > 0,x > 0), or copied from the last actively-updated 
row of cells if they are the physical "outflow" boundaries (the regimes y > 0, x < and 
y < 0,x > 0). For the purposes of current paper, the difference between these two BCs is 
not significant and resultant density profiles and torque calculation are almost identical (of 
the level of 10~ 3 ). We found that with the pure outflow BC fluid entering the simulation 
box accumulates some non-zero velocity perturbation on top of the pure linear shear velocity 
profile. This affects calculation of variables derived from the simulated velocity field, such 
as potential vorticity, which are analyzed in Paper II. We use the inflow/outflow y BC for 
our simulations. 



At the start of a simulation run we instantaneously turn on the potential of the planet in 
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the center of the linearly sheared fluid flow with uniform surface density. This gravitational 
perturbation immediately excites an inner and an outer density wave propagating away from 
the planet, as well as strong transients: vortices that appear near the planet but travel on 
horseshoe orbits away from it. Before measuring the wave properties, we run simulations 
for ~ 15 orbital periods to let these time-dependent structures move away from the planet. 
We also tried gradually turning on planetary gravity by linearly increasing streng t h of the 



poten tial from zero to its full value within several orbital periods (jMuto et al.ll2010l ; iLi et al. 



20091 ). but did not find this trick as effective at removing transient structures as simply 



waiting for them to move away from the planet. 



4. Result and comparison with linear theory 

We now present our main results and compare them with the predictions of linear 
theory. To mitigate possible nonlinear effects we consider only very low planetary masses 
in our calculations, ranging between 3.2 x 10~ 2 M th and 2.8 x 10~ 3 M th . According to Eq. 
(JT)) at 1 AU this mass range corresponds to M p between 0.4 M e and 3 Lunar masses (note 
that they are sound speed dependent). As the basis for comparison we primarily use the 
distribution of the perturbed density ( §4.11) and the evolution of the angular momentum flux 
carried by the density wave ( §4.2p . We carry out linear theory verification both in physical 
(coordinate) space and in Fourier space. 



4.1. The density wave profile 

The most basic and direct comparison between numerical simulations and theory can 
be performed using the spatial distribution of the perturbed fluid variables. In this work we 
choose the perturbed surface density ST, = T — T Q as the primary variable for comparison. 
Since ST is spatially concentrated along a narrow wake it makes sense to compare with 
theory both the overall shape of the wake in x — y coordinates, and the density distribution 
across the wake. 

We determine the wake shape in the following way. At each value of x we find the value 
of y, called y w (x), at which ST reaches its maximum value ST max (x). We then compare the 
run of y w (x) with the theoretical prediction ([3]), which applies in the shearing sheet geometry 
(GR01). The results are shown in Figure [2j The two curves agree well far from the planet, 
at |x| > h. The discrepancy at \x\ < h is not surprising because in this region (a) the density 
wave is not yet fully formed and (b) the dispersion relation ([I]) is significantly affected by the 



h 2 term, which causes t he wake profile to devia te from the theoretical prediction ([3]) valid 
far from the planet, see lOgilvie fc Lubowl (120021 ) for details. 



Next we investigate the evolution of the wave amplitude as it travels away from the 
planet by looking at the behavior of maximum amplitude 5T, max (x) = 5T l (x,y w (x)) as a 
function of x. We plot this quantity scaled by So(M p /M t h) (this normalization removes the 
dependence of 5T, max on E and M p ) in Figure [3^- 

At small x < h the value of <5£ max is large and decreases with x. This behavior has 
nothing to do with the density wave since the region in the immediate vicinity of the planet 
represents a quasi-static atmosphere that forms inside the planetary potential well right after 
the planetary gravity is switched on. This structure is rarely mentioned in the analytical 
calculations of planet-disk interaction but it shows prominently in realistic simulations. In 
hydrostatic balance (in the absence of background velocity caused by differential rotation) 
density profile in the vicinity of the planet should be given by (assuming isothermal equation 
of state with constant sound speed c s ) 

S atm (p) = Eoe-^WAS (18) 

In real disks, the background velocity field distorts the atmospheric profile from the circularly- 
symmetric (with respect to planet) shape predicted by this equation. This explains why 
Satm(^) shown Figure [3]3 slightly overestimates <5£ max at small x. 

The effect of the atmosphere on 5Tj max rapidly decreases beyond several h from the 
planet, and SH max starts to rise with increasing x. This is because density wave excitation 
occurs at |x| ~ h and from now on the behavior of 5T, max is determined by the wave-like 
density perturbation. Conservation of the angular momentum flux Fh carried by the wave 
forces the wave amplitude to increase with distance oc x 1 ^ 2 (GR01). As Figure shows this 
scaling agrees reasonably well with the numerically computed 5Y> max (x) far from the planet. 

Note that initially, at separations of several h from the planet, 5H max increases with x 
faster than the theoretical x 1 ^ 2 scaling, which is caused by the still ongoing accumulation of 
the angular momentum flux Fg at these separations. In other words, at this location Fg 
has not yet reached its asymptotic value given by Eq. see §4.21 for more details. 

Far from the planet, beyond x = Z s h (/ s h = 7-5h in the case shown in Figured), the 
peak <5£ rapidly goes down relative to the analytical x 1//2 scaling. This behavior is caused 
by the appearance of the shock at / s h and subsequent dissipation of the angular momentum 
flux carried by the wave. These nonlinear processes are explored in detail in Paper II. 

Finally, we examine evolution of the density profile across the wake as a function of 
radial separation x. We do this by making an azimuthal cut through the density field at 
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fixed x and shifting the resultant one- dimensional density profile in y by y w given in Eq. 
03]). To eliminate the increase of resulting simply from the angular momentum flux 
conservation we normalize the density perturbation by x 1 ^ 2 . We additionally scale 5E by 
E (M p /M t h), as in Figured- As mentioned before, in linear theory the normalized density 
profile should be independent of M p at a fixed separation x. This is indeed seen in Figure 
Ek, where we plot scaled profiles of computed for several values of M p at x — 1.33/t, and 
find no significant difference between them. The adopted value of x is small enough for the 
linear theory to apply even for the highest M p (rj 0.03M t h) used in making this Figure, i.e. 
x is always considerably smaller than (cf. Figure [SJo and the discussion in §4.41) . 

Previously, GR01 computed the evolution of the density profile as a function of x in 
the linear regime. In Figure [3b we provide an analogue of Figure 1 of GR01 by showing the 
scaled density cuts at the same values^ of x — 1.33/i, 2.67/i, Ah, 5.33/i — as in GR01. This 
particular calculation uses a very small planet mass M p = 3.7 x 10 -3 M t h corresponding to 
l s rs 8h, see Eq. Q, which puts the values of x used in making this Figure well inside the 
shock position and reduces the impact of nonlinearity on the density profile, see §4.41 

The effective width of the density profile in this Figure does not significantly vary with 
x and stays at the level of several h at all times. We note, however, that this property 
holds only for the density cuts passing through the wake at fixed x, as shown in Figure [SJo. 
On the contrary, the width of the density profile cutting through the wake at fixed y (not 
shown here) is shrinking as \x\~ l , which follows directly from the dispersion relation ([T]) since 
k r oc \x\ for x > h. 

Our numerical results agree with linear calculations of GR01 quite well, with quantita- 
tive differences at the level of 10% or less. In particular, we are able to reproduce several 
subtle features of the wake profile evolution found in GR01, which are highlighted by arrows 
in Figure [3b. These are (a) the increase of ST, with x in the density trough in front of the 
wake, (b) the decrease of <5£ with x in the density trough behind the wake, and (c) the slight 
drift of density profile with respect to theoretical wake position y = y w (x) given by Eq. ([3]) 
as x increases. There are also some minor differences between the Figure [3)d and the results 
presented in GR01, e.g. the lower (by about 10%) amplitude of far from the planet, at 
x > 3h, in our case (note that the vertical axis of [3b is different from the vertical axis in 
Figure 2 of GR01, since we use M t h and h instead of Mi and I = 2h/3 as mass and length 
units). We speculate that this may be explained by the gradually accumulating nonlinear 
effects in our simulations, which are absent in linear calculations of GR01. 



2 Note that GR01 normalized x by the "Mach l"distance I = (2/3)h rather than h as we do here. However, 
the physical values of x for the density cuts shown in Figure [3^ are the same as in GR01. 
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Previously, iMuto et al.l ( 120101 . Figure 6) have shown density profiles resulting from their 
2D hydro simulations in cylindrical geometry, which are similar to those presented in Figure 
[21 However, their calculations were carried out for planetary masses (minimum M p = 5 x 
10~ 2 M th ) larger than used in our work (ma ximum M p = 3_ L 2x 10~ 2 M th ), which does not allow 
them to isolate linear effects clearly. Also, IMuto et al.l (120101 ) did not attempt quantitative 
comparison of their numerical results with the linear theory predictions. 



4.2. The angular momentum flux and torque density 



We now investigate the behavior of the angular momentum flux (AMF) carried by the 
wave Fh(x) as a function of x. Linear calcu lations of this quanti ty in F ourier space have been 
first p erformed in GT80, and later refined in lArtymowiczl ( 1993al lbl) and lKorycansky fc Pollack 
(119931 ) . The torque density dTn I dx, wh ich i s a spatial derivative of Fh ( x) has been inferred 
from simulations in bate et al. I J2003h and b'Angelo fc Lubowi tooi . boioh . although 
direct comparison with the linear theory was provided in these studies. 



no 



The analytical expression ([6]) for the integrated angular momentum flux scales linearly 
with the planetary semi-major axis r p . In the shearing sheet geometry employed in our 
simulations r p is ill-defined, and it makes sense to redefine angular momentum flux as Fjj = 
Fjjjr v . Subsequently we will drop tilde for brevity and denote such normalized momentum 
flux as Fjj. We compute Fjj(x) from our data according to the following definition: 



00 

Fjj(x) = S J dy5v y v X} 



(19) 



where v x is the radial velocity of the fluid, and 5v y is the velocity perturbation with respect 
to the background shear profile in the azimuthal direction. 

We compute the torque density (normalized by r p ) exerted by the planet on the fluid 
at separation x as 



dT » - - [ 



dx J a dy ' 

—00 

and the integrated torque accumulated by the wave at x is just 

X 

T H (x) = / —^-dx. 
/ dx 



(20) 



(21) 



Conservation of angular momentum ensures that the derivative of Equation ( 1T91) (the slope 
of the AMF) coincides with the torque density (12"!]) in the absence of dissipation. 
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4-2.1. Comparison with theory in Fourier space 

We first compare our simulation results with linear theory in Fourier space by studying 
the behavior of the torque cutoff function FH{m)/F^ KB (m), see §2j Previously GT80 
computed variation of this quantity as a function of m, and the resul t is sho wn in Figures 



2 and 3 of their paper, providing a basis for comparison. lArtymowiczl ( Il993af ) subsequently 
refined this calculation by arguing that evaluating the planetary potential at separations 
slightly displaced from the classical Lin dblad resonance posi tion fl2]) should lead to more 



accurate results. Adopti ng this approach lArtymowicz ( 1993b ) provided a simple analytical 



prescription (Eq. (25) of lArtymowiczlll993bl ) to reproduce the behavior of Fh (m) / Fjj (m 



which can also be compared with our calculations. 

Substituting v x and 5v y in the form of their Fourier integrals in Eq. f fT9|) and ma- 
nipulating the resulting expression one can write Fh(x) as the integral over the azimuthal 
wavenumber k y : 



F H (x) = j dk y F H , k (x) } (22) 
o 

F H , k {x) = 47rS [Re (v X)k ) Re (Sv ytk ) + Im (v X)k ) Im (6v yjk )] , (23) 



where v Xjk and Sv y>k are the Fourier transforms of v x and 8v y , which themselves are functions 
of x. Even though the available analytical calculations of the torque cutoff function were 
performed only in the limit x ^> 1, we chose to retain the dependence of Fh(x) on x to 
explore the evolution of the harmonic content of the AMF with distance from the planet. 

Instead of the discrete theoretical WKB Fourier harmonics of the flux F^ KB (m) given 
by flS|) we use the following continuous version^: 

Fir = {k y hf (^e) 2 [2tfo(2/3) + *(2/3)] 2 , (24) 

where we also normalized the final expression by r p . 

We can now compute real and imaginary components of v Xyk and 5v Vjk using our nu- 
merical data and then obtain FH, k (x) from Eq. ( |23l) . Dividing the result by F^j^ B provides 
us with the ratio Fn,k{x)/ F^^ B , which can then be compared with the existing theoreti- 
cal calculations of the torque cutoff. This is done in Figure [5] where we plot the numerical 



3 Transition between the discrete and integral representations of the AMF is performed by replacing m 
with k v r p and summation over m with integration over r p dk y . 
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FH,k(x)/ '^Hk computed at different separations from the planet against the se mi-analytical 



calcula tions of the same quantity (in the limit \x\ 3> h) performed by GT80 and lArtymowicz 



( Il993bl ). These calculations were done for rather small M p = 1.2 x 10 2 M t h (corresponding 



to l sh w 5h) allowing us to see linear wave evolution out to large separations. 

One can see in Figure |5^l that the agreement between theory and simulations is generally 
quite good. Our resu lts seem to agree better with GT80 torque cutoff prescription than with 



Artymowica (jl993bl ). However, at the level of accuracy available to us we are not able to 
firmly discriminate between these two torque cutoff prescriptions. At intermediate values of 
\i = k y h « 0.5 — 5 our numerical results for Fn,k(x) / 'F^^ B are essentially independent of x 
and pass between the two aforementioned analytical torque cutoff prescriptions. Only weak 
sensitivity of F H>k (x)/ F™* 3 on x is expected since beyond \x\ = 2h the power in azimuthal 
harmonics corresponding to this range of k y h should have already been fully accumulated 
by the density wave, and no evolution should arise. Good agreement with the GT80 results 
is also illustrated in Figure [5b, which shows the torque cutoff function (multiplied by /i 2 ) in 
linear space, as in Figure 3 of GT80. 

The situation is different for smaller values of k y h < 0.5, where the numerical curves 
fall below the analytical asymptotic behavior Fu,k/ FjfJ^ B — > 1 as k y h — > 0. Also, there is 
noticeable evolution of Fn,k{x) / F^^ B with x, with nonzero power extending down to smaller 
and smaller values of k y as the wave travels further from the planet. This is because the 
low-order harmonics correspond ing to small k y h contribute t o the wave flux predominantly 



at separations larger than ~ h (IKorycansky fc Pollack! Il993l ). see Eq. (T5]). As a result, the 



further the wave travels, the more power gets collected by the wave from these lower-order 
azimuthal harmonics of the planetary potential, but the very small values of k y < 0.05 never 
contribute to the wave flux even at |x| — 5/i. 

Our numerical results also exhibit noticeable disagreement with theory and clear evolu- 
tion with x at large values of k y h > 5, and we comment on their origin in §4.41 



4-2.2. Comparison with theory in physical space 

We now study the behavior of the AMF and torque density in physical space. Under- 
standing the spatial distribution of the latter quantity is important for properly computing 
the AMF in disks with non- uniform distribution of S in the vicinity of the planet, and is 
thus important for understanding the early stages of gap formation by massive planets. 

Here we look at disks with uniform distribution of S. In Figure [6] we show the behavior 
of the AFM Fh(x) computed according to the definition flTTJT) and also of the accumulated 
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torque Th(x) given by Eqs. (|2T|) . We display our results for four different values of M p , 
normalizing F H (x) and T H (x) by S (GM p ) 2 /(/?,Cg) — according to linear theory the shape 
of the resultant curves should then be independent of M p . Theoretical scaling of the torque 
oc M 2 is indeed largely confirmed by this Figure, but see §4.41 for more details. 

To facilitate the comparison with the torque calculation, AMF curves have all been 
vertically shifted by a small offset to cancel the non-zero value of Fh{0) at x = resulting 
from the residual gas motion in the horseshoe region and a small intrinsic non-zero starting 
value of the AMF (Rafikov Sz Petrovich, in preparation). This brings them in perfect accord 
with Th (x) curves close to the planet, which is expected to be the case in the absence of 
dissipation.]. At large separations, beyond Z s h, the Fh(x) curve starts falling below the Th(x) 
curve, which is caused by shock formation and dissipation of the wave AMF past this point. 
Quite naturally, this effect is more pronounced for higher M p , e.g. M p = 3.2 x 10~ 3 M t h and 
M p = 1.2 x lCT 3 M t h, corresponding to smaller 1^. For smaller masses shocks occur outside 
the simulation box and essentially no difference between Fh{x) and Th{x) can be seen. The 
post-shock behavior of the AMF will be studied in Paper II. 

In Figure [7] we plot the torque density dT H (x)/dx (12"01 as a function of x. The overall 
shape of the curve is consistent with pr evious nume r ical c alcula tions of the same quantity 



in cy lindrical coordinates performed by iBate et al. I (120031 ) and iD'Angelo fc Lubowl (12008 



20101 ). As expected from linear theory, beyond w 2h torque density drops dramatically, and 
the AMF rises only weakly, with \x\, which is evident from the flattening of AMF curves 
in Figure O This justifies the localization of the wave excitation to the immediate vicinity 
(within ~ 2h) of the planet, as described in £j2j 

To the best of our knowledge, no quantitative analytical description of how Fh(x) should 
vary with distance from the planet exists in the literature, despite its potential importance for 
the gap opening problem. To fill this gap we first tried the following theoretical prescription 
Fj} R (x) motivated by the calculations of GT80: at each x we compute the lowest azimuthal 
wavenumber fM m i n (x) for which the location of the Lindblad resonance xi, m < x using Eq. 
(J2J) and then calculate the AMF as the sum of all Fourier contributions Fjj(m) given by Eq. 
fllj) that correspond to m > (r p /h)fi min (x). The result, which we denote by Fj^ R (x), is given 
by: 

oo 



4 This additionally confirms low levels of numerical viscosity in our runs. 
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. , 2h 

where we adopt the torque cutoff function F^(/i)/F^ KB (/x) calculated in GT80. Note that 
F^ R (x) is in fact independent of the adopted value of r p since its calculation can be rephrased 
fully in terms of fi rather than m. 

This prescription essentially assumes that the AMF contribution F^(m) corresponding 
to a particular azimuthal wavenumber m gets picked up by the density wave in a step- 
like fashion, solely at a single discrete location corresponding to the position of the m-th 
Lindblad resonance. Previously, this recipe was used in GT80 to compute the asymptotic 
behavior of the torque density in the limit \x\ > h, considering azimuthal harmonics with 
1 < m < r p /h. Here we e ssentially extend this prescrip tion to the case of arbitrary m and 



x\. Linear calculations by lKorycansky fc Pollack! ( 119931 ) show that this approximation is not 



very accurate, and AMF contribution due to each potential harmonic is in fact accumulated 
by the wave over an extended range of x. Thus we should not expect F^ R (x) given by 
Eq. ( J25l) to accurately represent the real coordinate dependence of the AMF in the linear 
regime. Nevertheless, this prescription still provides a useful reference point and we plot it in 
Figure E^a) (the dotted line, denoted by LR theory). We also plot in Figure [7] the theoretical 
prescription for the torque density dT}f^(x) / dx obtained by differentiating Eq. fl25|) with 
respect to x. 

One can see from these Figures that the simple-minded theoretical prescription fT^BT) 
overestimates the torque density at small x and underestimates it at larger x. While the 
areas under the numerical and theoretical curves in Figure [7J which represent the full accu- 
mulated AMF, are close to each other, the profiles of the two curves are quite different. As 
a result, F]^{x) initially rises faster than the numerical Fh(x) in Figure® but eventually 
the two asymptote to a similar final value. The small remaining difference is caused by both 
the intrinsic numerical inaccuracy and the fact that unlike the AMF theoretical curve the 
numerical AMF curve has been shifted to have a zero starting point. 

As the next level of approximation to the behavior of the Fh(x) and Th{x) we used the 
semi-analytical calculations of these quantities in the linear regime by Rafikov & Petrovich (in 
preparation). The run of corresponding Th(x) is displayed in Fig. [6] by the dash-dotted curve 
(the curve has also been shifted downward to cancel the non-zero starting value, and we use 
the label "Linear theory" to indicate Rafikov & Petrovich's work in all the figures) and clearly 
demonstrates good agreement between the simulation results and the semi-analytical linear 
theory. Theoretical torque density based on this T H (x) is shown in Fig. [7] and also agrees 
well with the numerical results, much better than the derivative of Fj^ R (x). This additionally 
emphasizes the point that assigning the planetary torque produced by a particular potential 
harmonic to a single location corresponding to the respective Lindblad resonance (as done 
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in Eq. fl25l) ) is not a very accurate procedure. 

Finally, we also compare our numerical results with another quantity that can be com- 
puted in the linear approximation, namely the "pseudo-AMF" fj introduced in GR01 as 



This quantity reduces to Fh{x) at large separations from the planet |x| > h, where ST 
is determined solely by the density wave, propagating away from the planet. However, 
for |x| < h the pseudo-AMF is strongly affected by the presence of the aforementioned 
atmosphere around the planet (see §4.ip . which is not a propagating density perturbation 
and cannot be associated with the density wave (this explains the prefix "pseudo-"). 

We compute the theoretical value fj(x) in the linear approximation using unpublished 
data from GR01 and also derive fj(x) from our numerical data using the definition (T2TI) . The 
comparison between the two is shown in Figure El for very small M p = 3.7 x l(T 3 M th . One 
can indeed see that at small \x\ the pseudo-AMF increases with decreasing \x\, analogous to 
ST in Figure [3^, which is explained by the presence of the atmosphere accumulated in the 
planetary potential well. Comparison with Figure [6] also shows that fj(x) indeed reduces to 
Fh(x) at large \x\. Most importantly, the agreement between the theoretical fj(x) derived 
from linear theory and the numerical fj(x) is very good in the whole range of x explored, 
which provides additional verification of linear theory of the planetary density wave evolution. 



In linear theory the dynamics of the density wave is independent of the adopted equation 
of state (EOS), and depends only on adiabatic sound speed of the gas c s . To check this 
property we have run several test simulations with two different EOS: isothermal EOS and 
EOS with 7 = 5/3, while keeping c s the same. The results do not show any significant 
difference between the runs with different EOS in linear stage thus confirming theoretical 
expectations. However, in the nonlinear stage a particular form of the EOS does become 
important and this is investigated further in Paper II. 



We now discuss the nonlinear effects that arise in our simulations even prior to the ap- 
pearance of the shock, i.e. at \x\ < Z s h- Nonlinear distortion of the wave profile is unavoidable 



oo 




(27) 



— oo 



4.3. Sensitivity to EOS 



4.4. Emergence of the nonlinear effects 
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even for the very low amplitude density waves and even during the mostly linear phase of 
their evolution. The rate of accumulation of the nonlinear effects (which eventually become 
strong and lead to shock formation) scales with planetary mass and is lower for low M p . 

This point is illustrated in Figure Hb, which is very similar to Figure and presents 
azimuthal density cuts through the wake for different values of M p , but now recorded further 
out from the planet, at x = 4h. While in Figure Hk, at x — 1.33/i, the density profiles for all 
M p are essentially overlapping and agree with the analytical profile from GR01, the situation 
at x = Ah is quite different. 

At this location the (normalized) wake profiles for the two smallest masses, M p = 
2.8 x l(T 3 M th (corresponding to l s « 8.9/t) and M p = 5.7 x 10 _3 M th (l s « Q.7h), still agree 
with each other and the semi-analytical linear calculation from GR01 quite well. However, 
already for M p = 1.2 x l(T 2 M th (l s » 5h) one can see a noticeable distortion of the profile 
compared to the linear solution, with the leading edge (at the left) becoming steeper and 
the profile peak shifting to the left, towards incoming fluid. The wave for even higher 
M p = 3.2 x 10 _2 M t h (l s ~ 3.4/t) has already shocked at x = Ah, which is clearly reflected 
in its density profile: its leading edge is essentially vertical, exhibiting a discontinuity in the 
fluid variables across the shock. Thus, the role of the nonlinear effects, which is measured 
at a given x both by the slope of the leading edge of the profile and the shift of its peak 
relative to the linear solution, progressively increases with the planet mass. 

Nonlinear effects also manifest themselves in a variety of other, more indirect ways. In 
particular, they explain the evolution of the torque cutoff function in Fourier space and its 
deviation from theoretical predictions at high k y h > 5 (see §4.2.11) . Figured clearly shows 
a growing amount of excess power at the high values of k y h as the wave propagates to larger 
x. This behavior is naturally accounted for by the nonlinear wave steepening, which causes 
the transfer of the AMF power in Fourier space from the low-fc^ to high-fc y harmonics. 

Nonlinearity also affects the AMF behavior in physical space. Close inspection of Figure 
[6] reveals that the numerical torque and AMF curves tend to asymptote to lower levels for 
higher M p , when nonlinearity is stronger. The difference in asymptotic values at large x 
between the largest and the smallest planet cases is about 5%. Variation of M p also affects 
the height of dTn{x) j 'dx curve at smaller x ~ h, increasing the peak value of dTn{x)/dx for 
lower M p , as shown in Figure [71 This explains the dependence of normalized AMF and Th{x) 
curves on M p in Figure [6j the area under the differential torque curve is slightly (by ~ 3%) 
higher for M p = 2.8 x l(T 3 M th (the think solid line) than for M p = 1.2 x l(T 2 M th (the thin 
solid line), explaining the difference seen between the asymptotic values of integrated torque 
Th(x) in Figures [6b and(6H. Spurious dissipation related to numerical viscosity cannot be 
blamed for this effect since its effect on the AMF should be independent of M p . 
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We suggest that this behavior may be explained (at least partly) by the stronger nonlin- 
ear evolution of the wave profile for higher M p . According to Eq. (120]) nonlinear distortion 
of away from the theoretically expected-value should have an effect (growing with M p ) 
on the calculation of the torque density dT H (x)/dx. Since for M p not too close to M th this 
profile distortion becomes significant only far from the planet, where the planetary potential 
is weak, the impact of this nonlinear torque modification on the AMF curves is not very dra- 
matic, but it is definitely non-zero. We also note that one expects simulations with smaller 
M p to exhibit better agreement with the theory since they are less affected by the nonlinear 
effects. However in reality this is not always the case because the optimal set of numerical 
parameters for achieving the best possible numerical result for different M p is different (see 



4.5. "Negative torque" phenomenon 



Close inspection of Figure [7| reveals an interesting feature in the behavior of the torque 
density at large x. The inset in this Figure clearly shows that dTn{x)/dx becomes negative 
beyond X- ~ 3.2h, which is at odds with linear theory since according to GT80 dTn(x)/dx 
should not change sign at any non-zero x. The contribution of this negative torque density 
to the total integrated torque T H is rather small (under one per cent level), but its existence 
presents a challenge to the results of GT80. Further investigation reveals that this "negative 
torque" phenomen on is present a lso in previous lower resol ution numerical calculations of 
Bate et al. I hooi . Figure 12) and b'Angelo & Lubowl (bood Figure 7). 



To understand the nature of this effect we varied the numerical parameters of our 
simulations. For a given M p = 1.2 x 10~ 2 M t h we tried varying softening length r s from h/16 
to h/32 at fixed resolution of 256/ h, and then varied resolution from h/64 to h/256 for a 
fixed r s = h/16. As Figure [7] shows this does not make the negative torque go away and 
dTu{x) / dx still changes sign at the same value of x_. Even more interestingly, changing 
the planetary mass M p to lower M p = 2.8 x 10 _3 M t h while keeping everything else the 
same (resolution of 256/ h and r s = h/32) also does not affect the position of x_, as the same 
Figure clearly shows. Analogously, experiments with different sizes of the simulation box and 
different prescriptions of gravitational potential show that the existence of negative torque 
and the position of x are independent of these parameters as well. Moreover, simulations of 
D'Angelo fc Lubowl (120081 . Figure 12) suggest that the negative torque phenomenon shows 
up only at low planetary masses (M p < 0.03 Mj w 10 M e in their case), which makes 
possible explanations based on nonlinear effects unlikely. 



This point is confirmed in Rafikov & Petrovich (in preparation) who demonstrate that 
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the sign change of the torque density at large separation from the planet is in fact a purely 
linear effect. It can be understood by carefully accounting for the complex spatial structure 
of the torque density produced by each azimuthal harmonic of the planetary potential, which 
goes beyond the simple calculation presented in GT80. 



5. Sensitivity of results to numerical parameters 

We now explore how the results presented in previous section and their agreement with 
analytical theory are affected by the variation of purely numerical parameters in our simu- 
lations. Table [T] lists the numerical parameters that we varied and the values we explored. 
Values corresponding to our standard case are indicated in boldface. 



5.1. Numerical solver and order of accuracy 

In our simulations we compare two different Riemann solvers — Roe's linearized solver 
(jRoelll98ll ) an d HLLC (Toro 1999) — with three different algorithms for the spatial recon- 



struction step ( Stone et al. 20081 ): second order with limiting in the characteristic variables 



(denoted 2 in this work), which is the predominant choice in literature in this field, third 
order with limiting in either the characteristic variables (3c), or in the primitive variables 
(3p). We find that Roe and HLLC solvers yield nearly identical results both in terms of the 
density profile and in terms of Fh(x) and Th{x) (differences are at the 0.1% level). 

The effect of using a different order of accuracy on the profile of the density wave is shown 
in Figure [9^, where we plot azimuthal density cuts obtained with different numerical settings 

Table 1. Parameter space of the simulations 



Parameters 3 



Range 



Riemann solver (flux function) used 

Order of accuracy 

Boundary conditions in y 

Resolution of the simulation (cells per h) 

Planetary potential 5 

Softening length 

Equation of states of the fluid 

Mass of the planet Mp/M^ 



Roe, HLLC 
2, 3c, 3p 
Outflow, Inflow/Outflow 

64, 128, 256 



,(2) a(4) $(6) 



p ! ! »p 

1/8, 1/16, 1/32 
7 = 1, 7 = 5/3 
3.2 x 1CT 2 , 1.2 x KT 2 , 5.7 x 10" 3 , 3.7 x 1CT 3 , 2.8 x 10' 



a See Section [5] for details. 
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at x = 1.33/i (together with semi-analytical profile computed by GR01 in the framework of 
linear theory) and x = Ah. One can see that density perturbation is rather insensitive to 
variation of the order of accuracy. 

The sensitivity of the AMF calculation to varying order of accuracy is illustrated in 
Figure [TOa. where the linear theory curve based on the semi-analytical calculation by Rafikov 
and Petrovich (in preparation) is overplotted for reference. One can see that the higher order 
of accuracy (3c and 3p) results in a slightly lower asymptotic value of torque and AMF, and 
the difference in the asymptotic value for the two between the 2 case and the 3p case is 
about 5%. Furthermore, 3c and 3p cases agree with each other on the position where the 
AMF curve starts to deviate from the torque, ~ Ah, which is still somewhat smaller than the 
theoretical prediction Z s h ~ 5h in this case due to the relatively low resolution. However, the 
second order of accuracy moves this point inward, which means that the angular momentum 
carried by the wave starts to dissipate earlier in this case. In particular, we find that using 
a second order accuracy solver compared with the third order accuracy solver has similar 
effect to decreasing the resolution by a factor of 2 (in term of advancing the displacement of 
the AMF-torque separation point, see §5.21 and Figure [TUb). 

Finally, we also note that calculations with 3p order of accuracy produce considerably 
noisier velocity field (this explains the noisy corresponding AMF curve in Figure [TOb ) than 
the other two cases explored. This difference is important for calculation of the velocity-based 
variables such as potential vorticity, which is used in Paper II as the means of shock detection. 
For these reasons in our current simulations we use Roe solver with 3c reconstruction. 



5.2. Resolution 

In Figure [9]d we show the effect of varying resolution of our simulations on the density 
profile. In general, increasing resolution from 6A/h to 256 jh improves the agreement with 
linear theory, but only slightly. Lower resolution simulations overestimate the amplitude of 
<5£ by just several per cent compared to 256//z simulations. Thus, for the study of the density 
wave profile our simulations essentially reach convergence in terms of resolution already at 
QA/h. 

In Figure [TUb we show that the AMF and torque calculations are more sensitive to 
resolution, especially when the nonlinear effects become important. Increasing resolution 
causes the asymptotic values of Fh(x) and Th(x) to decrease, and the difference in the 
asymptotic value for the two between the 256/ h case and the 6A/h case is about 10%. The 
key factor that clearly shows the downside of low resolution is the location, at which the 
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numerical AMF curve deviates from the numerical torque calculation. Calculations shown 
in Figure [10] use M p = 1.2 x 10~ 2 M t h, which according to Eq. (Q corresponds to shocking 
length Z S h ~ 5h. In our highest resolution simulations (256//i), the AMF curve starts to 
deviate from the accumulated torque curve precisely at l sh « 5h, as expected from theory. 
But as we reduce resolution, the radial separation at which Fh(x) starts departing from 
Th{x) moves closer to the planet violating the agreement with theory. In 64 /h simulations 
such departure starts already at 3.2/t, considerably closer to the planet than the nominal 
shock location Z s h- This behavior is caused by the higher level of numerical viscosity arising at 
low resolution, which leads to the dissipation of angular momentum carried by the wave prior 
to shock formation (we verified this point by conducting some runs with explicit viscosity in 
paper II, which show behavior qualitatively similar to our low resolution runs here). Thus, 
high resolution is crucial to correctly capture the details of the nonlinear wave evolution. 
We elaborate on this issue in Paper II. 



5.3. Planetary potential and softening length r s 

The rate at which a given smoothed potential converges to $x is important for the 
problem of density wave generation. Indeed, the amplitude and spatial distribution of fluid 
perturbation is in the end determined solely by the potential behavior. As a result, if $ 
is substantially different from $x i n the region where most of the torque is exerted by the 
planet, i.e. at |x| ~ 0.2 — 2h, see Figure [7J then one can expect significant spurious effects 
modifying the density wave properties. 

A specific form of the potential can have a two-fold effect on the properties of the density 
wave excited by the planet. First, there is a deviation of the potential from Newtonian, 
which directly affects the wave excitation at a given separation from the planet. But in 
addition, different potentials result in different structures of the quasi-static atmospheres 
(see §4.ip that get collected in the planetary potential well. This has an effect on the pressure 
distribution in the vicinity of the planet and can also affect wave excitation by modifying the 
dispersion relation and displacing the effective positions of Lindblad resonances. In practice 
disentangling these two effects based on the results of simulations may be non-trivial. 

In Figure [9fc we show that for a fixed r s = h/16 varying the order of the potential does 
not have a significant effect on the density wave profile. However, the AMF and torque 
calculations are more sensitive to the form of the potential, as Figure [TUb demonstrates. In 
particular, AMF in calculations with less accurate potential (e.g. $p ) is lower than it is in 
more accurate potential (5% for $ p 4 ^). We note that the seemingly better agreement with 

(2) 

theory in the $p case is an accidental phenomenon at this set of other numerical parameters. 
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For example, panel (d) shows that using softening length h/8 instead of h/16 could shift the 
curves downward by ~ 8%, so if switching to h/8 in panel (c) then the case would come 
to perfect agreement with the linear theory (see §5.41 for additional discussion). 

In Figure [9U we explore the effect of varying r s for a fixed form of the potential $ p 4 ^ . We 
discover that lowering r s results in a higher amplitude of the density perturbation. On the 
other hand, Figure [TUB shows that a lower r s leads to a higher numerical AMF and torque, 
and the difference in the asymptotic value for the two between the r s = h/8 case and the 
r s = h/32 case is about 16%. Again, the better agreement between the theory and the h/8 
case is accidental. Our calculations presented in §Huse $ p 4 ^ potential with r s = h/32. 



5.4. Summary of convergence study 



We conclude that the overall shape of the density wave profile is generally less sensitive to 
the variations of numerical parameters (unless these values are very different from the ranges 
explored in this work, see below) than the spatial dependence of the AMF or integrated 
torque. This emphasizes the importance of these integrated quantities as diagnostics of 
various subtle effects influencing the wave properties. 

In the course of this numerical exploration we have also discovered that the "optimal" 
value of a particular numerical parameter (e.g. softening length r s ) showing closest agreement 
with the linear theory depends on the values of other numerical parameters (e.g. resolution, 
type of the potential, etc.) and on M p . In other words, there is no universal best choice for 
each numerical parameter of simulations. In both Figure 191 and [TU| only the relative position 
of the curves representing different parameters matters, not their exact locations and the 
discrepancies between them and the linear theory. On the other hand, as we demonstrate in 
Paper II, high order of accuracy, high resolution, and highly accurate form of the potential 
with small softening length are critical to properly resolve the nonlinear stage of the wave 
evolution. 

Hydrodynamic simulations of the disk-planet interaction must often be run for many 
orbital periods. This is the case e.g. in studies of planetary migration or gap opening, in 
which substantial effects appear on timescales of hundreds to thousands of orbital periods 
(depending on M p ). This severely restricts the choice of resolution and softening length 
at which such simul ations can run. Typical resolutions in global disk simulations found in 



the literature (IKlev 



Li et al. 2009: Yu et al. 



1999 



Bryden et al 



2010 



1999 



Muto et al 



Nelson et al.ll2000l : iD'Angelo et al.ll200l 12003 



20101 ) range from 32 jh to l/h, and r s is usually 



taken to be between O.lh to h, in combination with a second order of accuracy solver and 
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the second (see Eq. ffl5l) ) or even lower order representation of the planetary potential. On 
the other hand, accurate description of the migration rate and gap opening may only be 
possible if properties of the density wave excited by the planet are properly captured by the 
simulation. 

To test how reliably the wave structure can be reproduced in global simulations we 
compare two sets of simulations in our local setting with M p = 1.2 x l(T 2 M th . We run one 
simulation with the typical numerical parameters for local shearing box simulations found 
in recent literature (typical global simulations in literature use even smaller resolution and 

(2) 

r s ): second order of accuracy, resolution 32/ h and $ p potential with r s = h/A. Another 
simulation is run with our fiducial parameters — third order of accuracy, resolution 256/ h, 
and $p 4 ^ potential with r s = h/lQ — and their results are compared in Figure [TT1 

There are noticeable differences between the two runs. The simulation with typical liter- 
ature parameters produces smoother density profile, which deviates from analytical solution 
by about 20%, see Figure [TTR At the same time, our fiducial run shows deviations from 
analytical profile only at the level of 1%. 

The difference is even more pronounced when we compare the AMF and torque behavior 
in Figure [TTb. Our fiducial simulation yields asymptotic values of Fh(x) and Th(x) which 
are within several per cent of the expected theoretical value (JBJ). As expected, Fh(x) agrees 
with Th(x) all the way until x ~ 5h, which is precisely the shock location for the M p used, 
according to Eq. 0. On the contrary, in the typical literature simulation AMF Fh(x) 
starts deviating from the integrated torque curve Th(x) already at x ~ 2.5h, significantly 
in advance of the expected shock position (a factor of 2). This indicates that dissipation 
and transfer of the angular momentum carried by the wave to the disk fluid start earlier 
in the typical literature case than they should in reality. The most likely explanation for 
this behavior is the high level of numerical viscosity in the typical literature run. Moreover, 
the absolute asymptotic values of Fh(x) and Tjg(x) disagree with the theoretical prediction 
for asymptotic Fh by ~ 23% (comparing with ~ 2% in our fidicual simulation), which is 
quite significant. We note that the differences between our fiducial and the typical literature 
simulations persist also in experiments with larger planet masses, M p < M t h, while the linear 
regime of wave excitation still holds. 

Underestimate of the wave angular momentum flux in the low-resolution simulations 
may result in an underestimate of the planetary migration rate in global simulations. In our 
shearing sheet setup we cannot investigate the effect of resolution on the planetary migration 
rate: by design one-sided torques exerted by the inner and outer parts of the disk on the 
planet exactly cancel each other, while the migration speed relies on the small imbalance 
between them. But the very fact that the one-sided torques in low-resolution case deviate by 
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tens of per cent from the high-resolution case is suggestive that the migration speed should be 
adversely affected by poor resolution at the same level. Furthermore, the relative imbalance 
between the one-sided torques can also be a function of resolution potentially exacerbating 
the problem. 

The discrepancy in the AMF should also affect the fidelity of gap opening by the planet 
in typical literature global simulations. Lower AMF carried by waves means that the planet 
is less effective at repelling gas away from its semi-major axis, which would require higher 
M p to open a gap. In ad dition, the gap opening process is sensitive to the spatial distribution 



of the AMF dissipation ( lRafikovll2002 bl ) . As a result, spurious dissipation of wave AMF prior 
to shock formation in low-resolution simulations and the accelerated AMF decay after the 
shock may introduce artificial effects in the gap opening picture. To summarize, any numer- 
ical studies of processes in which density wave production and dissipation plays important 
role need to use high order of accuracy, high resolution, and accurate representation of the 
planetary potential with small softening length. 



6. A numerical issue in planet-disk simulations 

In this section we describe a commonly ignored numerical problem that we discover in 
disk-planet interaction simulations. Namely, we find that if the time-step alt of the simulation 
is too large, the code cannot resolve the motion of the fluid in regions where the fluid 
experiences large gravitational acceleration. This issue applies to simulations in general, 
but it is especially problematic when the orbital advection algorithms is implemented (the 
FARGO algorithm, citealtmasOO). 



Massetl (120001 ) introduced a modification of the standard transport algorithm (Fast Ad- 
vection in Rotating Gaseous Objects, FARGO) for simulations of differentially rotating sys- 
tems, which significantly speeds up the calculations. By getting rid of the average azimuthal 
velocity when applying the Courant condition, this technique results in a much larger time 
step dt, which is limited by the small perturbed velocity, than the usual procedure, where dt 
is limited by the full fluid velocity dominated by the differential rotation. FARGO has been 



implemented in Athena by IStone fc Gardiner! (120101 ) . 



In a shearing box without planets, the dynamical time scale is fi -1 and is uniform 
throughout the box. However, when the planet is introduced, the dynamical timescale is 
spatially varying, and could be characterized by the free-fall time, which is the time that it 
takes a fluid element to fall on the planet assuming a constant acceleration at the grid point. 
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For gravitational potential in the form of Eq. ffT6l) . the free-fall time is: 



/ 2(p2 + r 2)5/ 2 

tff =U GM p (p* + 2.5rf) (28) 

Note that i// depends on a specific form of the gravitational potential, and a more smoothed 
potential results in a larger tff. The smallest free-fall time in the entire domain, which is 
the limiting timescale for the simulation, occurs at the grid points that are adjacent to the 
planet, which have the smallest p (lf/2 of the cell size r c in Athena; the separation where 
the smallest tff occurs also depends on the form of potential). In our simulations, we always 
keep r"l ^> r 2 c (usually r s = 8r c ), so the smallest free-fall time is: 

GM P < 29 > 

On the other hand, according to the Courant condition when varying the resolution the time 
step dt oc r c . So the ratio tff jS /dt scales with numerical parameters and M p as 

*//.* . ' s (30) 





To properly resolve fluid dynamics in the vicinity of the planet, the ratio tff >s /dt should be 
kept above a certain threshold. 

Our calculations without orbital advection and with dt determined by the Courant 
condition usually have tff )S /dt ~ 150. However when we turn on the orbital advection 
algorithm, dt typically increases by a factor of ~ 10, and the ratio tff jS /dt decreases by the 
same factor. We find that when the time-step is too large to properly resolve the fluid motion 
around the planet, the numerical results will be incorrect, as described below. 

Figure [T21 shows the azimuthal density profiles ST, at x = 1.33/i for a set of simulations 
using identical numerical parameters and orbital advection algorithm but with different 
tff. s /dt ratio, which we achieve by manually setting dt to be a fraction of the dt set by 
the Courant condition in FARGO. The two cases with highest tff_ s /dt yield the density 
profile in agreement with the theoretical prediction, also demonstrating the convergence of 
the result at high tff )S /dt. However, the density profiles in simulations with lower tff )S /dt 
clearly deviate from theory, with smaller tff yS /dt leading to larger discrepancy. In the run 
with the smallest used tff iS /dt = 18 which corresponds to dt set by the Courant condition 
in FARGO (representing maximum FARGO speed up), the density peak even becomes a 
density trough. By experimenting with Athena, we generally find that our simulations of 
the disk-planet interaction produce converged results agreeing with theory only when 

t fft ./dt> 100, (31) 
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while this critical number may depend on the numerical method used in a particular code or 
some other numerical conditions. We note that in principle all disk-planet simulations have 
to satisfy condition ( 1HT|) to ensure correct results, but in experiments we find simulations 
with the orbital advection algorithm tend to violate this condition much more easily than 
simulation without the orbital advection algorithm. 

Apart from incorrectly reproducing density perturbation, disk-planet simulations vio- 
lating the tff tS /dt constraint exhibit other serious problems. In particular, simulations with 
small tff s /dt do not achieve a steady state. We illustrate this in Figure IT21 by showing the 
density profile at half simulation time for two cases. While the tff, s /dt = 144 case reaches 
a steady state and develops time independent density profile at half time, the tff^/dt = 72 
case does not, and its density perturbation decreases with time. 

In simulations with tff a /dt lower then the critical value, the fluid element approaching 
the planet at small \x\, instead of moving on a horseshoe orbit, gets trapped by the planetary 
potential, leading to the formation of a fluid loop around the planet. As simulation progresses 
the rotational velocity of the fluid loop becomes higher and higher, and eventually the 
centrifugal force evacuates the region around the planet. This effect, which is a purely 
numerical artifact of the too small tff iS /dt, is visually very similar to (and can easily be 
confused with) the gap opening phenomenon, which is a real physical effect. Simulations 
exploring the gap opening process with large planet mass and low resolution are likely to have 
small tff s /dt violating condition (T5T1 see Eq. I5U1) . which may have detrimental effects on the 
results. We note that using Athena with orbital advection algorithm, numerical simulations 
with M p ~ M t h and a large r s = h/4 must have an effective resolution at least 64/ h in 
the vicinity of the planet to avoid this problem (a more smoothed potential may weaken 
this condition, through). However, again we emphasize that the critical value of tff tS /dt in 
simulations with the same numerical parameters but using different code and solver may be 
different from ours. 

The tff iS /dt threshold severely limits the ability of FARGO to speed up calculations 
of disk-satellite interaction. For most of our simulations, tff iS /dt is already rather close to 
the threshold value without orbital advection algorithm, so there is not much room left to 
increase dt, which is what the orbital advection algorithm aims to achieve (all our simula- 
tions presented outside this section are done without orbital advection technique). However, 
in many other studies of differentially rotating systems, such as the investigation of magne- 
torotational instability (MRI) in accretion disks, the point mass objects are absent and the 
characteristic dynamical time scale is always long enough for the orbital advection algorithm 
to be a extremely useful tool for speeding up the simulations. At last, we note that both dt 
and tff jS may vary during the simulation, and the constraint fl3Tl) on their ratio should be 
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satisfied in the entire domain throughout the simulation time. 

7. Summary 

We have conducted a series of hydrodynamical simulations to study the details of the 
gravitational interaction between low mass planets and a protoplanetary disk, and to test 
predictions of the linear theory of density wave evolution. Our simulations assume local 
shearing sheet geometry and are carried out in 2D at very high resolution to reduce the 
effect of numerical viscosity. We focus on both the excitation of the density waves and 
their propagation away from the planet in the linear regime. To mitigate the effects of 
nonlinearity we consider very small planetary masses, starting at 0.4 M ffi and going down to 
3 Lunar masses at 1 AU. 

We extract the spatial distribution of the density perturbation induced by the planet 
from our simulations and compare it with the predictions of linear theory. We find good 
agreement between the two, at the level of several per cent when high resolution (typically 
256 /h) is employed. We also investigate the spatial dependence of the angular momentum 
flux carried by the wave and the distribution of torque induced by the planet on the disk, 
again finding good agreement with theoretical predictions. In particular, we are able to 
reproduce the theoretical behavior of the torque cutoff in Fourier space. 

We also find various manifestations of nonlinear effects in our simulations even while 
the waves are formally linear. These include (1) the noticeable steepening of the density 
profile at larger values of M p far from the planet, (2) the slight variation of the asymptotic 
value of the AMF with planetary mass in addition to the expected F H oc dependence, 
causing the discrepancy with linear theory prediction at the level of several per cent, and (3) 
the growth of power in high azimuthal wavenumber harmonics of the AMF in Fourier space. 

By carefully studying the spatial distribution of the torque density dT H (x)/dx in our 
simulations we discover an interesting "negative torque" phenomenon: dT H (x)/dx changes 
sign at some radial separation from the planet (at \x-\ ~ 3.2h in our simulations), which 
contradicts the analytical results of GT80. This effect can however be understood in the 
framework of the linear theory as shown in Rafikov & Petrovich (in preparation). This 
feature of the torque distribution has only weak effect on the total accumulated torque in 
our simulations. 

We also carried out a detailed investigation of the effect of different numerical settings 
on our results in linear regime. We explored the influence on wave properties of (1) different 
Riemann solvers with different accuracy, (2) spatial resolution, (3) different forms of the 
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softened planetary potential, and (4) softening length of the planetary potential. We find 
the spatial distribution of the AMF and torque to be a more sensitive probe of the effects of 
various numerical parameters on the wave evolution than the distortions of the density wave 
profile. Based on this study we conclude that a very high resolution (ensuring low numerical 
viscosity), high order of accuracy, and an accurate prescription of planetary potential with 
small softening length are critical for accurately reproducing the key features of the wave 
evolution in linear regime. 

We demonstrate that low order of accuracy, low spatial resolution, and inaccurate po- 
tential with large softening length often employed in global disk-planet interaction studies 
can severely affect the fidelity of their results, especially in applications to planet migration 
and gap opening. Specifically, we conduct a test run with typical numerical parameters 
from recent literature. Comparing with the analytical theory, the test run produces a lower 
amplitude of the density wave (by ~ 20%), a lower final accumulated torque (by ~ 23%), 
and a largely advanced starting point of wave dissipation (a factor of ~ 2). 

Most of our own calculations were carried out with third order of accuracy, resolution of 
256//i, softening length r s = fo/32 and a potential that rapidly converges to Newtonian (as 
(r s /p) 4 ) at large separation p from the planet. This set of numerical parameters allows us to 
obtain excellent agreement with linear theory in Athena, but is not meant to be universal 
for other codes. However, the way in which we make the comparison with the theory and 
the agreement we achieve may serve as a standard framework for future code tests. 

We also discover a numerical problem which is largely ignored in previous simulations of 
disk-planet interactions. To follow the fluid motion correctly, the time-step in a simulation 
has to be smaller than the local dynamical timescale by a certain factor (~ 100 in our 
case using Athena) in the entire domain, including the region where the fluid experiences 
the largest gravitational acceleration (e.g. the vicinity of the planet). In the context of 
numerical disk-planet studies, violation of this condition leads to inaccurate calculation of 
wave properties, lack of convergence on long time scales, and spurious repulsion of gas 
by planet similar to gap opening. This timescale issue applies to disk-planet simulations 
in general, but it particularly affects the ones which use the orbital advection algorithm 
(FARGO), in which cases a significant increase of the time step is usually allowed to speed 
up calculations. 

In Paper II, we will continue our investigation of the disk-planet interaction by looking 
at the details of the density wave evolution in the nonlinear phase. 



31 



Acknowledgments 

We are grateful to Jeremy Goodman, Takayuki Muto, Zhaohuan Zhu, and an anonymous 
referee for useful discussions and comments. The financial support of this work is provided 
by NSF grant AST-0908269 and Sloan Fellowship awarded to RRR. 

REFERENCES 

Artymowicz, P. 1993a, ApJ, 419, 155 
Artymowicz, P. 1993b, ApJ, 419, 166 

Bate, M. R., Ogilvie, G. I., Lubow, S. H., & Pringle, J. E. 2002, MNRAS, 332, 575 

Bate, M. R., Lubow, S. H., Ogilvie, G. L, & Miller, K. A. 2003, MNRAS, 341, 213 

Bryden, G., Chen, X., Lin, D. N. G, Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 
344 

DAngelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 

DAngelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 

DAngelo, G., & Lubow, S. H. 2008, ApJ, 685, 560 

DAngelo, G., & Lubow, S. H. 2010, ApJ, 724, 730 

Dong, R., Rafikov, R. R., & Stone, J. M. 2011. larXiv:1109.2590l 

Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 

Gammie, C. F. 1996, ApJ, 457, 35 

Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509 
Gardiner, T. A., & Stone, J. M. 2008, Journal of Computational Physics, 227, 4123 
Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 
Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 (GT80) 
Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793 (GR01) 
Kley, W. 1999, MNRAS, 303, 696 



-32 - 

Korycansky, D. G., & Pollack, J. B. 1993, Icarus, 102, 150 

Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52 

Lin, D. N. C, & Papaloizou, J. C. B. 1993, Protostars and Planets III, 749 

Masset, F. 2000, A&AS, 141, 165 

Muto, T., Suzuki, T. K., & Inutsuka, S.-i. 2010, ApJ, 724, 448 

Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 

Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950 

Papaloizou, J. C. B., & Terquem, C. 2006, Reports on Progress in Physics, 69, 119 

Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., k Artymowicz, P. 2007, Proto- 
stars and Planets V, 655 

Rafikov, R. R. 2002a, ApJ, 569, 997 

Rafikov, R. R. 2002b, ApJ, 572, 566 

Rafikov, R. R. 2006, ApJ, 648, 666 

Roe, P. L. 1981, Journal of Computational Physics, 43, 357 

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 
137 

Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142 

Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 

Toro, E. F. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin: 
Springer) 

Ward, W. R. 1986, Icarus, 67, 164 

Ward, W. R. 1997, Icarus, 126, 261 

Yu, G, Li, H., Li, S., Lubow, S. H., & Lin, D. N. C. 2010, ApJ, 712, 198 



This preprint was prepared with the A AS IATgX macros v5.2. 




Fig. 1. — A typical image of our simulations, showing the density structure and the spiral 
waves. The quantity plotted here is S/E . M p = 2.09 x 10~ 2 M t h in this case. 
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Fig. 2. — Location of the peak density in the wake in the x — y plane measured from 
simulations (solid curve), compared to the analytical wake shape (dashed curve, Eq. (J3J). 
As expected, they agree far from the planet. Numerical peak density position is nearly 
independent of the simulation parameters. 
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Fig. 3. — Behavior of the density perturbation in a simulation with the following parame- 
ters: M p = 3.7 x lCT 3 M th (corresponding to ~ 3.6 Lunar masses at 1 AU), shocking distance 
l sh pa 8h. (a) Variation of the peak value of the relative density perturbation 5E/E with x 
[solid curve). Analytical profile ffl8|) of the quasi-static gaseous envelope collected inside the 
planetary potential well (without background shear) is shown by the dotted curve. Dashed 
line shows theoretical scaling <5£ oc x 1//2 (normalization is arbitrary) resulting from conser- 
vation of the angular momentum flux carried by the wave, (b) Azimuthal density profiles 
<5X (scaled by the planetary mass and normalized by (x/h) 1 ^ 2 ) at several values of x. To 
be compared with theoretical density profiles computed in linear approximation at the same 
locations in GR01 (their Figure 1). See text for details. 
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Fig. 4. — Azimuthal (normalized) profiles of the density perturbation for different M p (la- 
beled in panels) at two radial distances from the planet: (a) x = 1.33/i and (b) x = Ah. 
Analytical density profile from GR01 is shown by thick solid curve. See text for details. 



-37- 




Fig. 5. — Spectral decomposition of the angular momentum flux Fg(x) carried by the wave 
in azimuthal Fourier harmonics. Numerical results are based on isothermal simulation for 
M p = 1.2 x l(T 2 M th (corresponding to 0.14 M ffi at 1 AU and l sh rj 5h) with 128/h resolution, 
r s = h/16. (a) Ratio Fh ,k (x) /F^^ B (the so-called torque cutoff function) of the numerical 
Fourier spectrum of the AMF to the analytical WKB calculation of GT80 as a function of 
k y h, plotted at several values of x. For comparison we s how analogou s quant ity (in the limit 
\x\ —> oo) computed by GT80 (see their Figure 2) and lArtymowiczl (1l993al . their Eq. 25). 
The origin of excess power at high k y is discussed in §4.41 (b) The same quantity multiplied 
by fi 2 (using the x = 3h curve) and plotted on linear scale to facilitate comparison with 
Figure 3 of GT80. 
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Fig. 6. — Integrated torque Tjj(x) and the angular momentum flux carried by the wave Fjj(x) 
(shifted vertically so that Fh(0) = to simplify comparison with Th{x)) as a function of 
x. Different panels correspond to different values of M p (labeled on the plot, the same as 
in Figure H]). Analytical prescription 02 5)) for Tjj(x) motivated by the results of GT80 is 
shown in panel (a) by the dotted curve. The theoretical prediction based on the linear semi- 
analytical calculation (Rafikov and Petrovich, in preparation) is shown by the dash-dotted 
line in all panels. Both Fh{x) and Tjj(x) are scaled by S (GMp) 2 / (hc%), which should make 
their shapes independent of M p in the framework of linear theory. 
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Fig. 7. — Torque density dT H (x)/dx (Eq. [20]) as a function of x computed for two different 
values of M p and different numerical parameters. The thick dotted curve shows the deriva- 
tive of Eq. ( 1231) motivated by the results of GT80, and the thin dotted curve shows the 
linear theory result based on the semi-analytical calculation by Rafikov and Petrovich (in 
preparation). Torque density is scaled by S (G , M p ) 2 / (h 2 c^) removing the dependence on M p 
in linear approximation. For M p = 1.2 x 10~ 2 M t h we show results from three simulations 
with different adopted resolution or r s . An inset zooms in on a region near x = 3h, where 
we find the torque density changing sign. The location (x_ = 3.2h, shown by arrow) at 
which this happens is found to be insensitive to variations of simulation parameters or M p 
(M p = 2.8 x 10~ 3 Mth was also explored). 
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Fig. 8. — "Pseudo-AMF" defined by Eq. (l27|) from our numerical simulations (solid curve) 
and semi-analytical linear calculations of GR01 (dashed curve). The agreement between the 
two is very good in the linear regime (for M p = 3.7 x 10~ 3 M t h used in simulation wave shocks 
at Z s h ~ 8h). 
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Fig. 9. — Results of the numerical convergence study. Different panels show results of 
simulations (normalized azimuthal density cuts, analogous to those shown in Figure HJ at 
x = 1.33/i and x = 4h) in which a single numerical parameter was varied: (a) order of 
accuracy for the solver, (b) resolution, (c) softened planetary potential, (d) softening length 
of the potential. The standard set of numerical parameters is M p = 1.2 x 10~ 2 M t h, $p 4 ^ 
potential, r s = h/lQ, 128//i resolution, and 3c accuracy. Results of the simulation with 
this set are shown by solid line in all panels. The semi-analytical density profile from linear 
calculations of GR01 is shown only at x = 1.33/i to guide the eye. 
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Fig. 10. — Same as Figure El but for the spatial behavior of the (properly scaled to remove the 
dependence on M p ) integrated torque Th(x) (the thin line of each line type) and AMF carried 
by the wave Fh(x) (the thick line of each line type; shifted vertically so that Fh{0) = 0). 
The linear theory curve based on the semi-analytical calculation by Rafikov and Petrovich 
(in preparation) is shown for reference. 
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Fig. 11. — Comparison between a simulation with our fiducial numerical parameters (solid 
line, third order of accuracy, resolution 256 /h, $p potential with r s = h/lQ), and a simula- 
tion with numerical parameters typically adopted in recent planet-disk simulations (dashed 
line, second order of accuracy, resolution 32/ h, potential with r s = h/A). We use 
M p = 1.2 x 10~ 2 M t h (l s h ~ 5/i) for both runs, (a) Azimuthal (scaled) density cut (similar to 
that in Figure [9]) at x — 1.33/i showing lower amplitude of the density perturbation in the 
typical literature run. (b) AMF and integrated torque as a function of x (similar to Figure 
ITUT) for both runs. The amplitudes of the density perturbation, F H (x), and T H (x) are lower, 
and the AMF dissipation starts earlier in the typical literature case, compared to our fiducial 
simulation and the linear theory. The linear theory curve is based on the semi-analytical 
calculation by Rafikov and Petrovich (in preparation). 
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Fig. 12. — Azimuthal density profiles (scaled by the planetary mass and normalized by 
(x/h) 1 / 2 ) at x — 1.33/i, for simulations using orbital advection algorithm (FARGO) with 
different tff }S /dt ratio (we manually set dt in different cases to be a fraction of the value of 
dt set by the Courant condition in FARGO). Dashed lines show the density profiles at half 
simulation time for two representative runs to illustrate whether temporal convergence has 
been achieved. Theoretical prediction (GR01) is also plotted to guide the eye. Simulations 
are done with resolution 64//i, potential with r s = h/8, and M p = 3.2 x 10 _2 M t h. 



